


# This scripts loads the intervention-country level data, fits the ideal point model, and saves the stan fit object and extracted country ideal points
# Inputs: (1) ar5ar6intervention_country.csv
# Outputs: (1) iso_names.RData
#           (2) stan.fit.RData
#          (3) idealpts.csv
# Note: This script takes approx 7 hours to run on a Macbook pro with 6 cores



# rm(list = setdiff(ls(),"root")) #clear workspace
# Set Working Directory
path <- 'replication_final'
setwd(path)

# Load Libraries
library(tidyverse) # version ‘2.0.0’
library(rstan) # version ‘2.32.6’
library(countrycode) # version ‘1.6.0’
rstan_options(auto_write = TRUE)



#--------- Stan Code
stan.code <- "
    data {
        int<lower=2> K; // Number of response categories
        int<lower=0> N; // Number of observations
        int<lower=1> D; // Dimensionality of ideal points
        int<lower=1> I; //Num Countries
        int<lower=1> J; //Num Interventions
        int<lower=1, upper=I> ii[N]; //Country for observation n
        int<lower=1, upper=J> jj[N]; //Intervention for observation n
        int<lower=1, upper=K> y[N];
        int<lower=1> Hidx; // Index of country with positive ideal point
        int<lower=1> Lodx; // Index of country with negative ideal point
    }
    parameters {
        row_vector[D] beta[J]; // discrimination parameter of item j
        row_vector[D] x[I]; // ideal point of country i
        ordered[K - 1] c[J]; // item-specific cutpoints
    }

    model {
        for (d in 1:D) {
            to_vector(x[,d]) ~ normal(0, 1);
            x[Hidx,d] ~ normal(1, 0.0001); // Prior for Hidx country with positive ideal point
            to_vector(beta[,d]) ~ normal(0, 5);
        }
        x[Lodx,1] ~ normal(-1, 0.0001); // Prior for Lodx country with negative ideal point
        for (m in 1:J) {
            c[m] ~ normal(0, 25);
        }


        // Logit 
        for (n in 1:N) {
            y[n] ~ ordered_logistic(dot_product(x[ii[n]], beta[jj[n]]) , c[jj[n]]);
        }
    }"




#--------- Load Data
# Intervention-Iso3c Level Data
statements<- read_csv(paste0(path,'/ar5ar6intervention_country.csv'))
statements

# Restrict to inrterventions with atleast 2 countries
subsetsinterventions<- statements %>% 
                        group_by(interventionid) %>%
                        summarise(n = sum(abs(x), na.rm = T))  %>% 
                        filter(n>=2) %>%
                        ungroup()  %>% 
                        pull(interventionid)

# create roll call dataframe
rcdf <- statements %>% 
                filter(interventionid %in% subsetsinterventions) %>% 
            arrange(statementid, interventionid, iso3c) %>%
            dplyr::select(statementid, iso3c,interventionid, x) %>% 
            group_by(iso3c) %>%
            mutate(iso3c_id = cur_group_id())  %>% 
            ungroup()  %>% 
            group_by(interventionid) %>%
            mutate(interventionid_id = cur_group_id())  %>%
            ungroup()  %>%
            group_by(statementid) %>%
            mutate(statementid_id = cur_group_id())  %>%
            ungroup() 


K <- 3 # Response Categories
D <- 1 # Dimensionality of Ideal Points
II <- rcdf %>% pull(iso3c)  %>% unique() %>% length() # Number of Countries
JJ <- rcdf %>% pull(interventionid) %>% unique() %>% length() # Number of Interventions
SS <- rcdf %>% pull(statementid) %>% unique() %>% length() # Number of Statements
N <- II*JJ # Number of Observations
y <- (rcdf$x)+2 # Response Variable
iso_names<- rcdf %>% distinct(iso3c) %>% pull(iso3c)

# save vector of (ordered) iso3c names
save(iso_names, file = paste0(path, "/iso_names.RData"))

i_n <- rcdf$iso3c_id # Country for observation n
j_n <- rcdf$interventionid_id # Intervention for observation n
s_n <- rcdf$statementid_id # Statement for observation n
Hidx<- which(iso_names=="DEU") # Set Germany as high anchor
Lodx<- which(iso_names=="SAU") # Set Saudi Arabia as low anchor


# put in a stan data list
stan.data <- list(K=K, N=N, D=D, 
                  I=II, J=JJ, 
                  ii=i_n,jj=j_n, 
                    y=y,
                  Hidx=Hidx, Lodx = Lodx)


#--- Fit Stan Model 
stan.fit <- stan(model_code=stan.code, 
                data=stan.data, 
                iter=2000, warmup=1000, 
                chains=4, thin=1, 
                verbose=TRUE, cores=6, seed=1234)

save(stan.fit, file=paste0(path, "/stan.fit.RData"))

# load(file = paste0(path, "/processing/stan.fit.RData"))
#---- Get summary of Stan model
stan.out<- summary(stan.fit)$summary



#----- Extract Ideal Points
x.est<- bind_cols(
    x1up = stan.out[str_detect(rownames(stan.out),"x\\[\\d+,1]"),'97.5%'],
    x1down = stan.out[str_detect(rownames(stan.out),"x\\[\\d+,1]"),'2.5%'],
    x1=stan.out[str_detect(rownames(stan.out),"x\\[\\d+,1]"),'mean'],
    iso3c = iso_names
) %>% 
    arrange(x1) %>%
    mutate(order1 = 1:nrow(.)) 

x.est$country <- countrycode(x.est$iso3c, "iso3c", "country.name")
x.est$country[x.est$iso3c=="EUU"] <- "European Union"
x.est<- x.est %>% mutate(iso3c = ifelse(iso3c=="EUU", "EU", iso3c))
x.est$iso3c <- factor(x.est$iso3c, levels = x.est$iso3c)



write_csv(x.est, paste0(path, "/idealpts.csv")) # Save ideal points




